home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / ode-initval / cstd.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  4.5 KB  |  185 lines
  1. /* ode-initval/cstd.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. #include <config.h>
  21. #include <stdlib.h>
  22. #include <math.h>
  23. #include <gsl/gsl_errno.h>
  24. #include <gsl/gsl_math.h>
  25. #include <gsl/gsl_odeiv.h>
  26.  
  27. typedef struct
  28. {
  29.   double eps_abs;
  30.   double eps_rel;
  31.   double a_y;
  32.   double a_dydt;
  33. }
  34. std_control_state_t;
  35.  
  36. static void *
  37. std_control_alloc (void)
  38. {
  39.   std_control_state_t * s = 
  40.     (std_control_state_t *) malloc (sizeof(std_control_state_t));
  41.   
  42.   if (s == 0)
  43.     {
  44.       GSL_ERROR_NULL ("failed to allocate space for std_control_state", 
  45.                       GSL_ENOMEM);
  46.     }
  47.  
  48.   return s;
  49. }
  50.  
  51. static int
  52. std_control_init (void * vstate, 
  53.                   double eps_abs, double eps_rel, double a_y, double a_dydt)
  54. {
  55.   std_control_state_t * s = (std_control_state_t *) vstate;
  56.   
  57.   if (eps_abs < 0)
  58.     {
  59.       GSL_ERROR ("eps_abs is negative", GSL_EINVAL);
  60.     }
  61.   else if (eps_rel < 0)
  62.     {
  63.       GSL_ERROR ("eps_rel is negative", GSL_EINVAL);
  64.     }
  65.   else if (a_y < 0)
  66.     {
  67.       GSL_ERROR ("a_y is negative", GSL_EINVAL);
  68.     }
  69.   else if (a_dydt < 0)
  70.     {
  71.       GSL_ERROR ("a_dydt is negative", GSL_EINVAL);
  72.     }
  73.   
  74.   s->eps_rel = eps_rel;
  75.   s->eps_abs = eps_abs;
  76.   s->a_y = a_y;
  77.   s->a_dydt = a_dydt;
  78.  
  79.   return GSL_SUCCESS;
  80. }
  81.  
  82. static int
  83. std_control_hadjust(void * vstate, size_t dim, unsigned int ord, const double y[], const double yerr[], const double yp[], double * h)
  84. {
  85.   std_control_state_t *state = (std_control_state_t *) vstate;
  86.  
  87.   const double eps_abs = state->eps_abs;
  88.   const double eps_rel = state->eps_rel;
  89.   const double a_y     = state->a_y;
  90.   const double a_dydt  = state->a_dydt;
  91.  
  92.   const double S = 0.9;
  93.   const double h_old = *h;
  94.  
  95.   double rmax = DBL_MIN;
  96.   size_t i;
  97.  
  98.   for(i=0; i<dim; i++) {
  99.     const double D0 = 
  100.       eps_rel * (a_y * fabs(y[i]) + a_dydt * fabs(h_old * yp[i])) + eps_abs;
  101.     const double r  = fabs(yerr[i]) / fabs(D0);
  102.     rmax = GSL_MAX_DBL(r, rmax);
  103.   }
  104.  
  105.   if(rmax > 1.1) {
  106.     /* decrease step, no more than factor of 5, but a fraction S more
  107.        than scaling suggests (for better accuracy) */
  108.     double r =  S / pow(rmax, 1.0/ord);
  109.     
  110.     if (r < 0.2)
  111.       r = 0.2;
  112.  
  113.     *h = r * h_old;
  114.  
  115.     return GSL_ODEIV_HADJ_DEC;
  116.   }
  117.   else if(rmax < 0.5) {
  118.     /* increase step, no more than factor of 5 */
  119.     double r = S / pow(rmax, 1.0/(ord+1.0));
  120.  
  121.     if (r > 5.0)
  122.       r = 5.0;
  123.  
  124.     if (r < 1.0)  /* don't allow any decrease caused by S<1 */
  125.       r = 1.0;
  126.         
  127.     *h = r * h_old;
  128.  
  129.     return GSL_ODEIV_HADJ_INC;
  130.   }
  131.   else {
  132.     /* no change */
  133.     return GSL_ODEIV_HADJ_NIL;
  134.   }
  135. }
  136.  
  137.  
  138. static void
  139. std_control_free (void * vstate)
  140. {
  141.   std_control_state_t *state = (std_control_state_t *) vstate;
  142.   free (state);
  143. }
  144.  
  145. static const gsl_odeiv_control_type std_control_type =
  146. {"standard",            /* name */
  147.  &std_control_alloc,
  148.  &std_control_init,
  149.  &std_control_hadjust,
  150.  &std_control_free};
  151.  
  152. const gsl_odeiv_control_type *gsl_odeiv_control_standard = &std_control_type;
  153.  
  154.  
  155. gsl_odeiv_control *
  156. gsl_odeiv_control_standard_new(double eps_abs, double eps_rel,
  157.                                double a_y, double a_dydt)
  158. {
  159.   gsl_odeiv_control * c = 
  160.     gsl_odeiv_control_alloc (gsl_odeiv_control_standard);
  161.   
  162.   int status = gsl_odeiv_control_init (c, eps_abs, eps_rel, a_y, a_dydt);
  163.  
  164.   if (status != GSL_SUCCESS)
  165.     {
  166.       gsl_odeiv_control_free (c);
  167.       GSL_ERROR_NULL ("error trying to initialize control", status);
  168.     }
  169.   
  170.   return c;
  171. }
  172.  
  173. gsl_odeiv_control *
  174. gsl_odeiv_control_y_new(double eps_abs, double eps_rel)
  175. {
  176.   return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 1.0, 0.0);
  177. }
  178.  
  179.  
  180. gsl_odeiv_control *
  181. gsl_odeiv_control_yp_new(double eps_abs, double eps_rel)
  182. {
  183.   return gsl_odeiv_control_standard_new(eps_abs, eps_rel, 0.0, 1.0);
  184. }
  185.